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Abstract 

Aims. We discuss the applicability and reliability of the shapelet technique for scientific image analysis. 

Methods. We quantify the effects of non-orthogonality of sampled shapelet basis functions and misestimation of shapelet parameters. We 
perform the shapelet decomposition on artificial galaxy images with underlying shapelet models and galaxy images from the GOODS survey, 
comparing the publicly available IDL implementation with our new C++ implementation. 

Results. Non-orthogonality of the sampled basis functions and misestimation of the shapelet parameters can cause substantial misinterpre- 
tation of the physical properties of the decomposed objects. Additional constraints, image preprocessing and enhanced precision have to be 
incorporated in order to achieve reliable decomposition results. 
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1. Introduction 

Image analysis is the crucial technique and the prime objective 
in all observational sciences. Computer-based image analysis is 
supposed to provide us with information that human perception 
is not or only hardly capable of extracting, due to either the vast 
amount of data or the required accuracy. The decomposition of 
imaged objects into an orthogonal function set and analysis of 
the expansion coefficients is the preferred way of conducting 
image analysis since sufficient computer power is available. 



Refregier d2003l) proposed the 'shapelets' function set, 
composed of a scalable version of the eigenfunctions of 
the harmonic oscillator in quantum mechanics. They form 
an orthonormal basis set of Gauss-Hermite polynomials. 
Because of their continuity, finite extent and their smooth 
peaks, they offer themselves for decomposing galaxy im- 
ages or the like. In particular, they were proposed as an im- 
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galax ies dKellv & McKavl I2004I) and Sun spots (lYoung et al 



2005 ), and also in the field of medical computer tomography 
dWeissman etai]|2004l) . 



So far, the shapelet technique has proven to reconstruct the 
decomposed objects from the set of expansion coefficients in 
a way which look s visually good ; i.e. differences to the origi- 
nals are marginal ( Refregierl2003 ). Massev & Refregiei (2005) 



defined a goodness-of-fit measure to quantify, what a good re- 
construction means: the residuals have to be at noise level. For 
achieving that, the decomposition has to be optimized with re- 
spect to a set of parameters. As the parameters space is highly 
degenerate, additional constraints are necessary to choose a 
particular point in this space. But it is yet unknown, if the selec- 
tion of these constraints or the accuracy with which their fulfill- 
ment is verified introduces an uncertainty in or even a bias on 
the expansion coefficients. In fact, every study of object prop- 
erties derived from shapelet coefficients might be affected by a 
yet unclear error contribution. 

Another object of concern is the computational complexity. 
Since the shapelet method is known to be slow compared to 
other image analysis techniques, it is important to find ways to 
speed up the execution. 

This paper is organized as follows: In sect. 2, we summa- 
rize the basic relations for the shapelet function set and describe 
the procedure for finding optimal decomposition parameters. 
In sect. 3, we discuss potential problems that can arise from 
the optimized decomposition procedure. In sect. 4, we show 
how these problems can be remedied by means of additional 
constraints and image preprocessing. In sect. 5, we compare 
the design choices, the decomposition results (of artificial and 
observed galaxy images), the errors made by and the computa- 
tional performance of two shapelet implementations, the pub- 
licly available IDL code and our independently developed C++ 
code. We conclude in sect. 6. 
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2. Shapelet basics 

Following the work by lRefrepeil(l2003l) . the shapelet decompo- 
sition allows us to approximate a two-dimensional object (for 
example a galaxy image whose brightness distribution is given 
by I(x), centered at x c ) by a finite series 



"i 1 "z " max 

/(x) - / rec0 (x) = ^ InB a (x - x c ;B), 



(1) 



where x = (xi,Xz) and n = (ni,n 2 ). The two-dimensional 
shapelet basis functions 



B n (x;J3)=JT 1 (f> ni (jr 1 x 1 ) (f>„ 2 (B- l x 2 ), 



(2) 



are related to the one-dimensional Gauss-Hermite polynomials 

<t> n (x) = [2 n 7rinlT? H n (x)e-^, (3) 

with H„(x) being the Hermite polynomial of order n. 

These definitions imply that the shapelets' operating range 
is limited by 



: B(n max + 1)5 and 6 min = B(n max +1) '-, 



(4) 



denoting the maximal and minimal size of features in the im- 
ages that can be faithfully described by a decomposition using 
B and n max . 

After the decomposition, one can derive all information 
that depends on the brightness distribution of an object in the 
much smaller shapelet space instead of the real space, thus sav- 
ing memory and computation time. The equation relating the 
shapelet coefficients to the fl ux (or other quan tities such as the 
centroid position) is given by lRefre gier (2003), 
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Appropriate relati ons fo r the q uadrupole moments Q/j have 
been published bv lBerga d2005l) . 
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from w hich we can derive a com plex ellipticity measure fol- 
lowing (Bartel mann & Schneiderl2001 ). 



Gn - 222 + 2iQ n 
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(7) 



See (Ref regierl 120031) for further details about the shapelet 
method. 



Optimized decomposition procedure 

Equation ([TJ shows that a shapelet decomposition depends 
on four external parameters: the scale size B, the maximum 
shapelet order n max and the two components of the centroid 
position x c . The essential task for achieving a faithful shapelet 
decomposition is finding optimal values for these parameters, 
such that the residual between the original image and its recon- 
struc tion from shapelet coeffici ents is minimized. 
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Mass ev & Refregierl d2005l) defined a goodness-of-fit func 
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x c ) is a pixel vector of 



where R(J3, n max , x c ) = / 
the residual between the actual image brightness / and its re- 
construction by a shapelet model I reco , and V is the covariance 
matrix of the pixels. In the case of Gaussian noise with standard 
deviation cr„, V = cr^l. 

The number of coefficients is related to n max via Eq. (Q~|): 



n. 



coeffs — 
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(9) 



X 1 is normalized to the number of degrees of freedom and is 
equal to unity when the residual reaches the noise level. In this 
case, the decomposition procedure was able to extract all sig- 
nificant physical information present in the image. 

Since Eq. (JHJ is linear in the unknown shapelet coefficients 
/„, we ca n solve analyt ically for their values when^- 2 is mini- 
mal (e.g. lFriedenll 19831) . 



In 



{M T V- l My l M T V~ l I, 



(10) 



where the matrix M = M, ; (y6, n max , x c ) gives the value of the ith 
shapelet basis function sampled at pixel j, and / is a pixelized 
version of the brightness distribution I{x). 

Thus, optimizing the decomposition means finding the set 
of parameters for which x 1 becomes unity. One has to con- 
sider, though, that n max is a discrete parameter, which forbids 
using minimization algorithms for continuous parameters, but 
in turn restricts the parameter space severely. In addition, one 
must investigate whether the parameter set can be determined 
uni quely. 



Massev & Refregierl (120051) suggested the following proce- 
dure: Starting with n max = 2, the value of B is searched where 



dB\ 



0, 



(ID 



using a one-dimensional simplex minimizeiQ. From the 
shapelet coefficients, the deviation of the shapelet model from 
the given x c is computed, which is then subtracted from x c such 
that this deviation disappears. Then, with constant B, n max is in- 
creased until x 1 becomes unity or flattens (i.e. x 2 changes less 
with B than the uncertainty in x 2 )- At this new n max , an opti- 
mal value for B is searched with the simplex minimizer, again 
readjusting x c after each iteration. n max is then reset to 2 and 



1 See for example IPress et af. I <2002h for a description of the 
method. 
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increased again, using the already optimized values of f3 and 
x c , to possibly find a point in the parameter space with lower 
rimax- If this succeeds, (3 is optimized further. The procedure 
terminates if either n max or (3 converges. 

3. Potential shapelet problems 

X 2 as defined in Eq. ([8]l is a slowly varying continuous func- 
tion for the majority of images (see Fig.[4]i. This allows an easy 
and fast decomposition of objects, yielding reconstructions that 
look visually good and have low x 2 - The question arises if the 
decomposition is in all cases able to extract the physical infor- 
mation correctly. We show two classes of potential problems, 
where physical properties of the decomposed objects are in- 
deed incorrectly measured. One class is related to the violation 
of the orthonormality condition of the shapelet basis functions, 
the other one regards misestimation of the shapelet parameters. 



3.1. Orthonormality 

The matrix M T M is the covariance matrix of the expansion co- 
efficients, which in case of a complete orthonormal basis set 
is equ al to the identity matrix. As pointed out by iBerrv et all 
(2004J), the covariance matrix can well differ from the identity 
matrix and can lose the diagonal form or even its full rank. 
This happens just because a set of discrete basis vectors - de- 
rived from a continuous complete orthonormal basis set by 
sampling the continuous functions at points of a (regular) grid 
- can lose orthonormality, orthogonality or even completeness 
depending on the definition of the grid. Loss of orthogonal- 
ity introduces covariances among the expansion coefficients; 
losing completeness eliminates the possibility of inverting the 
covariance matrix in Eq. ( TfOl i, so that the x 1 approach is not 
applicable anymore. 



Adopting the methodology of IBerrv et al.l ( 120041) . we ex- 
plore the domains of non-orthonormality by inspecting the di- 
agonal and non-diagonal entries of the covariance matrix. 

3.1.1. Undersampling 

In Fig.Q] the largest and the smallest diagonal elements of the 
covariance matrix diverge at ft < 2 because the cell-size of 
the grid (1 pixel) is too large to represent the variations of the 
continuous shapelet basis function with n max = 10 and /? < 2. 
As can be seen in Fig. [2] the basis set also loses orthogonality 
in that domain, since the largest non-diagonal element of the 
covariance matrix is of the same order as the diagonal elements. 

To make things worse , undersampling is essentially equiva- 
lent to random sampling dBerrv et al.l l2~004). Since oscillations 
of the basis functions appear on smaller scales than the grid 
spacing, small shifts of the grid points can lead to arbitrary dif- 
fer ences in the function v alues. 

Massey & Refregieri J2005b suggested a way of dealing 



with undersampling: instead of using vectors sampled at cer- 
tain grid points, the value from integrating the basis functions 
within each pixel should be used. But this approach also leads 
to non-orthogonal basis vectors, independent of the scale size 
(see dotted curve in Fig. |2j. 




Figure 1. Largest (dashed, blue) and smallest (solid 
agonal elements of the covariance matrix M T M as 
of ft at Umax = 10. Deviation from unity indicates 
orthonormality. The image size was 50 x 50 pixels 
(25,25). 
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Figure 2. Largest non-diagonal element of the covariance ma- 
trix M T M as a function of ft at n„ Mf = 10 (solid, black). 
Deviation from zero indicates the loss of orthogonality. The 
second curve (dotted, red) shows the effect of using integrated 
basis functions. The third curve (dashed, blue) shows the effect 
of including a constant function in the basis set. The image size 
was 50 x 50 pixels with x c at (25, 25). 



3.1.2. Boundary effects 

Another problem arises when the scale size is too large for 
the shapelet basis functions to be contained inside the image. 
Since the shapelets need infinite support for their orthonormal- 
ity, power will be lost due to truncation at the image bound- 
aries, as shown by the lower curve in Fig. [1] at high /3. Again, 
also the orthogonality is violated in this domain (see. Fig. [2]). 




P. Melchior et al.: Reliable 




Figure 3. Example galaxy from GOODS CDF South. The im- 
age was chosen because of its typical deep field signal-noise- 
ratio and its significant substructure. The image size is 64x64 
pixels. 

3.1.3. Non-orthogonal elements 

Massev & Refregierl d2005l) mention the possibility of fitting 
the sky background brightness by adding a constant function 
to the set of the shapelet basis functions. As can be seen from 
the dashed curve in Fig. [2] this again violates orthogonality 
globally. This means, it introduces covariances between the co- 
efficients even in such domains of p where the shapelet basis 
functions themselves remain orthogonal. 
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ft max 

Figure 4. log 10 (y 2 ) for the decomposition of the galaxy from 
Fig. [3] The centroid was fixed. The dashed line delimits the 
X 1 < 1 region, where the shapelet parameters are degenerate. 

components of the ellipticity originates from the orientation of 
the galaxy along the top-right to bottom-left direction. 

Obviously, setting the maximum order to an arbitrary value 
will cause a misestimation of the coefficients and the derived 
quantities if the selected n max was too low to represent the entire 
object or too high to remain unbiased by surrounding noise. 



3.2. External parameters 

As mentioned in sect. [2] the set of external shapelet parame- 
ters will not be uniquely defined (see the degeneracy region 
with^ 2 < 1 in Fig . HI) . We need to specify which set should be 
chosen. This choice will affect the shapelet coefficients and the 
quantities derived from them. Quantifying the impact of param- 
eter misestimation on the results of the shapelet decomposition 
is the aim of this section. 

For this purpose, a visually selected galaxy from the 
GOODS CDF South (Fig. was decomposed until x 2 wa s 
compatible with unity (with minimal n„ mx ). Then, starting from 
the optimal values {n opt ax = 8,p op < = 5.39, x c °p' fixed from the 
image; cf. Fig. [4j, the decomposition was repeated with varied 
external parameters. For each decomposition, the flux and the 
ellipticity (defined by Eq. (|7]i) were derived from the shapelet 
coefficients together with the x 2 of the fit. 

3.2.1 . Variation of n max 

The maximum order n max was varied between and 20. In each 
case, the centroid was fixed, and p was chosen to minimize x 2 
at the given n max according to Eq. (fTTb . As expected, the fit 
improves with increasing n max (top panel of Fig . [5j ■ 

As can be seen in the middle panel of Fig. [5] the flux will 
be underestimated for low n max due to the lack of substruc- 
tures represented in the reconstruction. On the other hand, the 
reconstruction tends to pick up smaller noise features farther 
away from the center when n max exceeds the preferred value 
n °max- Thus the flux and especially the ellipticity (bottom panel 
of Fig. [5J become noisy at high n max . The correlation of the two 



3.2.2. Variation of/3 

In order to see how the accuracy of the minimizer (for finding 
the optimal P) or an arbitrary choice of the scale size influences 
the coefficients, the scale size was varied within ft 01 " — 2 < 
P < p°P ! + 2, corresponding to a change of * 40%. The other 
parameters were kept fixed. 

From the top panel of Fig. |6j one can conclude that x 2 is 
not strongly affected by a change in /?: The largest values are 
X 2 ~ 1.03. The shapelet decomposition is obviously able to 
cope with a mispredicted scale size and yet provides visually 
good-looking reconstructions. Minimizing x 2 w.r.t. p is fortu- 
nately simple due to the lack of local minima and the smooth- 
ness of the curve, but crucial because of the non-trivial depen- 
dence of the other quantities on p. 

The misestimation of the flux (middle panel of Fig. [6]l for 
P + P" pt is easily understood: Since the central peak is most 
significant, the peak height is essentially fixed for each recon- 
struction. If p < p opl , the reconstruction peaks sharply, falls off 
too fast and misses the outer parts of the object, thus the flux is 
underestimated. If p > p opl , the central peak becomes broader 
and the outer regions of the reconstructions are too bright, thus 
the flux is overestimated. 

The variation of the ellipticity estimator (bottom panel of 
Fig. |6j is also problematic. Since increasing p at constant n max 
increases 9 max (cf. Eq. (0]l), we see a similar behavior as before 
in the bottom panel of Fig. [5] where n max was increased: Due 
to the intrinsic orientation of the galaxy, both components of 
the ellipticity are correlated. With p < p opl , the model is more 
compact and the ellipticity is dominated by the galactic core 
which has excess flux top-right of the center. With p > P" 1 ", the 
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Figure 5. Impact of the variation of n max on the decomposition. 
In the top panel the x 1 of the fit is plotted. The other panels 
show the deviations of the flux (middle), and the components of 
the ellipticity (bottom) from the values at the selected optimum. 
The vertical line indicates the optimum value. 



outer parts become more important; because of excess flux far 
below the center the e, flip sign but remain correlated. 

3.2.3. Variation of x c 

In order to clarify if the determination of the centroid can safely 
be done during the iterations of the optimized decomposition 
without biasing the outcome, both components x c were varied 
by 5 pixels, thus moving the centroid along a straight line from 
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Figure 6. Impact of the variation of {! on the decomposition. 
All panels as explained in Fig. [5] 



(jc j - 5, x 2 ~ 5) to (x j + 5, x j + 5). The other parameters 
were kept fixed. 

As the top panel of Fig.|7]shows, the goodness-of-fit is not 
affected very much; largest values are x 2 ~ 1 -04. Interestingly, 
the x 1 minimum is not at the optimal value (determined by 
calculating the centroid in real space), but rather shifted by 
roughly one pixel towards the top right corner of the im- 
age. This indicate s that the optimized decomposition proce- 
dure (suggested by lMassey & Refregierl t 20051) and outlined in 
sect. I2J tends to converge to a slightly different centroid. In this 
case, the fit improves (lower x )> but it no longer represents 
the imaged object because the centroid is readjusted under the 
assumption that the object (without noise) can be perfectly rep- 



6 



P. Melchior et al.: Reliable Shapelet Image Analysis 




Figure 7. Impact of the variation of x c on the decomposition. 
All panels as explained in Fig. [3] 



resented by the employed shapelet model. This assumption is 
not satisfied in general. 

The underestimation of the flux is < 8% (middle panel of 
Fig. 13). If the center is misaligned, the flux is underestimated 
due the loss of regions in the reconstruction far away from the 
new center. The rise of AF at high Ax c indicates a positive flux 
region on the top-right side of the galaxy, which can be con- 
firmed by inspecting Fig. [3] 

The ellipticity estimator is highly problematic and rather 
unpredictable (bottom panel of Fig. [7]), because its deviation 
depends on the orientation of the object w.r.t. the image bor- 
ders. 



4. Solutions 

In this section, we show that none of the problems mentioned 
above is fundamental. They can all be remedied by introduc- 
ing additional constraints or image preprocessing. The general 
idea is that a more careful way of dealing with the shapelet 
decomposition, knowing its subtleties, will yield more reliable 
results. 

The first three solutions tackle problems related to the non- 
orthogonality of the basis functions, the last three show how 
optimal shapelet parameters can be found and how a realistic 
error estimate for the shapelet coefficients can be performed. 

4.1. Singular Value Decomposition 

Berr y et al. d2004 suggested to use a Singular Value 
Decomposition (SVD) of the matrix M to derive the shapelet 
coefficients: 



M = UI.V T => /„ = Vl T U T I, 



(12) 



where, if M is an mxn matrix, U and V are orthogonal matrices 
with the dimensions m x m and n x n, respectively. S is an m x n 
matrix whose diagonal elements are the non-negative singular 
values; all other entries are zero. £ contains only the inverses 
of the non-zero diagonal elements of S. 

This approach accounts for the non-orthogonality of the ba- 
sis set, caused by e.g. severe truncation at the image boundary 
or inclusion of a non-orthogonal element. The drawback of the 
SVD approach is its memory consumption: Since M is a matrix 
with dimensions n p i xe i s x n coe ff s , the matrix U has the dimen- 
sions Kpixeis x n pixels- For typically sized images this amounts 
to matrices sized of the order of 100 MB, rendering the imple- 
mentation inacceptably slow for most purposes. 

4.2. Noise subtraction & image segmentation 

Most images in astronomy and other fields are contaminated by 
noise. If the noise has zero mean, the shapelet decomposition 
can readily be used. If not, the additional noise background 
should be subtracted before the decomposition instead of being 
fit by a constant function. SExtractor or similar software can 
be used to estimate the statistical noise measures for the image 
(they will be needed for calculating x 1 anyway) and subtract 
the noise mean from the image brightness. 

Since potentially more than one object is contained in the 
image, the overall image has to be segmented in a next step 
such that each frame contains only one object. SExtractor or 
equivalents accomplish this by grouping pixels above certain 
brightness thresholds. If other objects overlap with this frame, 
they have to be masked with noise whose features were deter- 
mined in the previous step. Otherwise the locality condition for 
the shapelet decomposition - meaning that the decomposed ob- 
ject is centered at x c and has a limited extent - will be violated 
and the physical information obtained from this 'combined' ob- 
ject will be corrupted by the appearance of the other objects. 

Note that no noise-removal step has to (and should) be done 
by arbitrarily discriminating between the object of interest and 
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the surrounding noise. The great benefit of the shapelet tech- 
nique is that the optimized decomposition itself will provide 
this discrimination better than any simple object-finding algo- 
rithm, since it reproduces the actual shape of the object instead 
of e.g. convolving the image with predefined kernels. In fact, 
procedures similar to the shapelet decomposition (e.g. wavelet 
decomposition) have be en propose d to tackle this problem in a 
more quantitative way (Bertin l2005l) . 

4.3. Geometrical constraints 

Although the x 2 method for finding the shapelet coefficients 
(Eq. (fTOl i) works also in the case of a non-orthogonal set of 
basis vectors (then introducing covariances among the coeffi- 
cients), the cleanest way of dealing with the problems of un- 
dersampling and boundary truncation is to avoid them. 

Very faint and small objects that are decomposed with very 
low n„ mx and (5 < 1 must be rejected. Fortunately, these objects 
carry little physical information. 

The decomposition of objects with high S IN must use large 
nmax- Since the oscillations of the basis functions can then ap- 
pear on sub- pixel scales, their sampli ng becomes essentially 
random. As lMassev & Refregierl ( 12005 ) suggested, one can get 
rid of this by applying the additional constraint 



(13) 



meaning that the oscillation 'wavelength' should be larger than 
the grid spacings (1 pixel). If f3 is unrestricted, this effectively 
limits Umax via Eq. (0|. 

The opposite case of Q max becoming large compared to the 
image dimensions can easily be prevented by placing the ob- 
ject inside a frame that is large enough. This can simply be 
done during the image segmentation step. The additional bor- 
der size should be proportional to the object's extent, have a 
lower bound of the order of 10 pixels and be made such that the 
frame becomes square, since the extent of the shapelets is iden- 
tical along both dimensions. If the object is close to the image 
boundary such that the frame border cannot be created from the 
image pixels alone, the missing pixels should be set according 
to the local noise characteristics. The inclusion of the additional 
frame border also helps to remedy the degree-of-freedom sin- 
gularity: The form of Eq. © demands the additional constraint 



Wpixels tlcoeffsi 

which again limits n„ wx via Eq. (|9). 
4.4. Centroid independence 



(14) 



As shown in sect. 13.2.31 leaving x c as a free parameter of the 
decomposition can lead to deviations of the preferred shapelet 
centroid from the centroid of the object. Since the reconstruc- 
tion should represent the object as closely as possible, the cen- 
troid should be fixed to the value calculated from the object 
in real (pixel) space, instead of the one derived from shapelet 
coefficients, which will only yield an approximation of the cen- 
troid under the assumption that the true galaxy morphology is 



correctly represented by the employed shapelet model. In gen- 
eral this assumption does not hold, e.g. galaxy profiles often 
have steeper cusps than the best fit shapelet model. 

Apart from the fact that fixing the centroid does not assume 
a particular galaxy model, it eases the minimization ofx 2 w.r.t. 
P because any recentering induces distortions the simplex algo- 
rithm has to compensate in order to converge. 

We will show in sect. 15. 3 1 that the error of the centroid de- 
terminati on with our approach and the iterative recentering ap- 
proach of iMassev & Refregierl (12005b is « 0.5 pixels (we em- 
ployed shapelet models for the galaxies there) which is the ex- 
pected uncertainty due to pixelization of the image. 

4.5. Convergence of the optimizer 

As described in sect. 13.2.11 the minimization should not be 
stopped before it reaches a residual at noise level 



x 2 <i 



min(n max ) 



(15) 



with n max being as small as possible (for the sake of uniqueness 
of the parameter set). 

This can be easily realized by not applying any limit on 
n max . Unfortunately, the two constraints (Eqs. ( fT3l & ( TBI ) im- 
pose such a limit in special cases, leaving y 2 > 1 ■ The recon- 
struction can well be visually good, but not all physical infor- 
mation can be extracted or trusted. These objects should there- 
fore be classified as 'silver', in contrast to the 'golden' ones 
which obey Eq. dl~5b . For a follow-up analysis, one could then 
decide whether the 'silver' sample should be used or the selec- 
tion should be restricted to the 'gold' sample. 

Equation (T5[ also states that n max has to be minimal. It is 
thus insufficient to find a point in parameter space with^ 2 * 1, 
but one has to ensure that this parameter set is also the one with 
the lowest possible n max . In practice this requires additional it- 
erations at lower n max . Condition (fT~5T > thus also prefers Eq. (fTTT i 
to be fulfilled, but this has to be ensured separately. 

A flattening condition for limiting n max when x 2 does not 
improve significantly should not be used because it would lead 
to an incomplete optimization when^ 2 is not strictly monoton- 
ically decreasing with n max , e.g. due to a faint satellite at the 
object's boundary. 

4.6. Extended error estimation 

So far, the only source of errors mentioned to give rise to coef- 
ficient errors was the pixel noise (via V in Eq. (fTUli). If the basis 
set remains orthogonal, there will be no additional uncertainty 
coming from the x 2 method. 

Additional - fortunately independent - errors are intro- 
duced by the lack of knowledge on the optimal parameter set. If 
the centroid position is fixed, one can assume its error to vanish 
(which is not exactly true, see sect. 15.31 ). The same is true for 
the error of n mclx when the minimizer converged, obeying con- 
dition dl5t . The simplex minimizer for finding ft (fulfilling Eq. 
( fill ) stops when it has localized f3 within a previously defined 
interval. Therefore, f3 will only be known up to an uncertainty 
A/3 that should, of course, be very small, maybe on the order 
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of 1%. Additional uncertainty in fi is introduced by pixel noise 
and pixelization, so that we expect the actual A/3 to be some- 
what larger. 

The question then arises how A/3 can be translated into 
an uncertainty in the shapelet coefficients. The answer to this 
question can be obtain ed by using the rescaling operation, 
given in appendix A of (Refregier 2003), which describes the 
change of the shapelet coefficients due to a change of the scale 
size. Assuming fi will for sure be within a 10% interval (this 
assumption will be justified in sect. 15.3b , one could identify the 
change of each individual coefficient by corresponding rescal- 
ing with a 3-cr error. In contrast to the error introduced by un- 
correlated pixel noise, the uncertainty in fi will in general affect 
coefficients differently. As we will also show in sect. 15.31 this 
error will be the dominant contribution to the coefficient error. 



5. Implementation comparison 

All of the above mentioned weaknesses may become impor- 
tant, but it is not a priori clear in which cases they may play a 
decisive role. We thus proceed from a description of the prob- 
lems to a comparison of the decomposition results of the pub- 
licly available shapelet implementation^!, written in the inter- 
preted programming language IDL, and our independent im- 
plementation in C++. We also compare errors made during the 
decomposition and the performance of the two codes. 



5.1. Design choices 

The publicly available implementation of the shapelet tech- 
nique, which has already been used for a considerable num- 
ber of studies (see sect. [TJ, is written in IDL, which implies 
several drawbacks: it is not efficient in dealing with large nu- 
merical problems, and licenses are quite expensive. Thus, we 
decided to reimplement the shapelet decomposition indepen- 
dently, in C++ using Open Source Software only. Since C++ 
is not full-featured for numerical computations - in contrast to 
IDL - our implementation links with very powerful external 
libraries: GNU Scientific LibrarjQ, ATLASQ lapack| and 
Boosj§ The dependence on external libraries has the minor dis- 
advantage that they have to be compiled and installed before 
our code can run, but takes full advantage of the development 
and improvement effort spent on these libraries. Using C++ 
and the numerical libraries, our code is roughly one order of 
magnitude faster than the IDL implementation (see sect. I5.4I 
for the performance comparison). 

Furthermore, we decided to define our implementation as 
a C++ library, such that other codes that need the shapelet 
method can easily link with it. This approach is not possible 
with IDL because IDL scripts are interpreted and therefore not 
executables in their own right. 



In particular, our implementation performs the following 
steps to decompose images, employing the solutions pointed 
out in sect. [4] For the task of estimating the noise character- 
istics and segmenting the image into frames, SExtractor is 
the standard choice in astronomy. Unfortunately, it is not avail- 
able as library, only as standalone executable. Since we did 
not want our library to use system calls, we decided to im- 
plement algorithms similar to those in SExtractor in our li- 
brary. Thus, the whole image preprocessing (noise estimation 
and subtraction, image segmentation and removal of overlap- 
ping objects) is done internally in our code, before the object's 
frame is passed to the shapelet decomposition. Also, the cen- 
troid of the object is directly computed from the cleaned frame. 

For the determination of the shapelet coefficients, we use 
the x 1 method (see sect. [2]); the SVD method is also imple- 
mented, but not used as default. We optimize the parameters fi 
and Umax, which have to fulfill the constraints given by Eqs. ( fTTT i 
and ( fT3l > - (fTBT) . We pay particular attention to find the lowest 
possible n max . x c remains fixed at the value derived from the 
cleaned frame. After the convergence of the minimizer, the er- 
rors on the shapelet coefficients are computed from the pixel 
noise and the uncertainty in /3 (see sect. 14.6b . If the decompo- 
sition violates condition (TT3T >. the object is classified as 'silver' 
so that it can be excluded from the further analysis. A flattening 
condition is not employed. 

5.2. Decomposition results 

We now investigate the first decomposition results from 1,000 
simulated images with underlying shapelet models, and then 
the results derived from the decomposition of 2,660 galaxies 
from the GOODS survey. 

In both cases, in order to guarantee the comparability of 
the implementations, the images were cut to contain the object 
of interest in a square frame, large enough not to run into the 
problem of boundary truncation - the IDL code was modified 
such that it does not further cut the frame. The noise character- 
istics were estimated using our noise algorithm and then passed 
to both shapelet codes. We can thus focus on the optimization 
procedure employed by the different codes. 



5.2.1. Simulated images 



2 http : //www . astro . caltech . edu/~r jm/shapelets/ 
We used the most recent stable version 2. 1/3. 

3 http : //www . gnu . org/ sof tware/gsl/ 

4 http : //math- atlas . sourcef orge . net 

5 http://www.netlib.org/lapack/ 

6 http : //www .boost . org | 



Kelly & McKay! ((2004) identified the ten most powerful (i.e. 
largest on average) shapelet coefficients of galaxy images from 
SDSS. Using these coefficients and their variances, we de- 
fined a multivariate Gaussian probability distribution. To create 
more individual galaxy shapes, we included minor coefficients, 
whose variance was chosen to be smaller than of the major 
ones; their mean was set to zero. Sampling from this proba- 
bility distribution, we generated 1,000 flux-normalized galaxy 
images with n max = 8 and 2 < /3 < 10. To each of these shapelet 
models, we applied a moderate level of Gaussian noise with 
zero mean and constant standard deviation of 2 • 10 . 

We then decomposed these images into shapelets using the 
IDL code and our C++ implementation with n max < 20 and 
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Figure 8. Comparison of the final decomposition x 2 for tne 
two implementations using simulated images with underlying 
shapelet models. 

1 < < 20. These restrictions are quite loose, but helpful to 
constrain the parameter space for this comparison. 

Generally speaking, both implementations are able to de- 
compose the given images with a x 2 near unity as required 
(see Fig. [HJ. The IDL version typically needs higher orders 
to achieve the same goodness-of-fit: The mean n max for C++ 
is 9.79, but 11.97 for IDL (see Fig.|9]l. Both implementations 
tend to higher n max than the underlying shapelet model because 
noise can obscure the smallest features of the model most ef- 
fectively, requiring higher additional shapelet modes to correct 
for that. 

The other optimized shapelet parameter is the scale size /3. 
In Fig. [10] the scale size of the decomposition results is com- 
pared to the one of the underlying shapelet model. For both 
implementations one can clearly see a correlation, with signif- 
icantly larger scatter in the IDL case. Interestingly, the images 
for which the IDL decomposition misestimates the scale size 
considerably form the group of images with n max = 20, the 
maximum allowed order. This shows that a misprediction of 
the scale size is compensated by using higher shapelet orders. 

To quantify the overall quality of the decomposition, we 
computed the mean and the Pearson correlation coefficient for 
several quantities of the model with the decomposition output 
(see Table [TJ. It becomes obvious from these numbers that the 
decomposition results are more reliable in the C++ than in 
the IDL implementation, regarding both the value of the ex- 
ternal shapelet parameters and the quantities derived from the 
shapelet coefficients. 

5.2.2. GOODS images 

We selected 2,660 galaxies from GOODS subject to the con- 
straint that the object frame does not contain more than 10,000 
pixels. The reason for this selection is twofold: The speed of 
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Figure 9. Comparison of the final decomposition n max for the 
two implementations using simulated images with underlying 
shapelet models (n max = 8). Since n max is a discrete parameter, 
the data points have been randomized with a Gaussian distribu- 
tion (standard deviation 1/3 pixel). 
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Figure 10. Comparison of the final decomposition /?' ++ (bot- 
tom, black) p' dI (top, red) with /3° of the underlying shapelet 
models. The values of /3 ldl have been shifted upwards by 5 units 
for better discrimination. 



the code scales at least linearly with the number of pixels, so 
we wanted to keep the procedure reasonably fast, and the typ- 
ical objects for shapelet decomposition (e.g. in weak lensing 
studies) are faint and rather small. 

The noise mean and variance were estimated using our 
noise algorithm. The noise mean was subtracted from the im- 
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Figure 11. Comparison of the final decomposition^ 2 for the Figure 13. Comparison of the final decomposition y3 with 
two implementations using GOODS galaxy images. fi° ++ for the selection of GOODS galaxy images. 
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Figure 12. Comparison of the final decomposition n max for the 
two implementations using GOODS galaxy images. The num- 
bers have again been randomized with a Gaussian distribution 
(standard deviation 1/3 pixel). 



ages before decomposition; the noise variance was passed to 
the shapelet codes. 

As can be seen from Fig.QTl both implementations are able 
to achieve a reconstruction with x 1 ~ 1 ■ There is a clear cutoff 
at^ 2 + cr(x 2 ) for the C++ code that is less prominent for the 
IDL code. This could be due to the use of a flattening condition 
during the optimization procedure in the IDL code, which al- 
lows higher^ 2 if the residuals do not reduce significantly when 
increasing n max . 



The comparison of the maximal orders shows a clear corre- 
lation which has to be expected from the two implementations, 
since it shows that both use low order for small and faint ob- 
jects and higher order for brighter and bigger objects (see Fig. 
rT~2l> . Again, with IDL there is a larger fraction of objects at the 
nmax = 20 limit, but the trend is not as clear as with the simu- 
lated images. 

The most surprising comparison is the one of the final scale 
sizes (see Fig. [131 . We can see the expected spread, but there 
is an alignment of data points along constant fir 1 . For some 
reason, the IDL optimization procedure prefers some values of 
fi, although this is a free and continuous parameter. This effect 
could explain the generally larger scatter of fi' dl for the simu- 
lated images, too. As we have discussed in sect. 13.21 the result 
of such an inconsistency will affect the whole optimization out- 
come and therefore the quality and reliability of any follow-up 
analysis based on shapelet coefficients. 

5.3. Decomposition errors 

With the simulated galaxy images and the decomposition re- 
sults of the two codes at hand (cf. section [5".2.1t , we can inves- 
tigate which effect gives rise to what kind of error. 

It seems reasonable to assume that the uncertainty in the 
shapelet parameters percolates through the decomposition pro- 
cess and creates the scatter in the quantities derived from 
shapelet coefficients. But the misestimation of the shapelet pa- 
rameters can also create systematic biases: Figs.[6]&[7]clearly 
show that it is quite likely to underestimate the flux F if fi is un- 
derestimated or the centroid x c is shifted away from its optimal 
value; exactly this behavior is conspicuous in the IDL results. 

When we compare the distribution of centroid determina- 
tion errors AR C = |Ax c | = |x c - x°| (see Fig. [14] and last line 
of Table Q]), we can clearly see that both codes are affected by 
a mean uncertainty of w 0.5 pixels, which is compatible with 
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Table 1. Means and Pearson correlation coefficient r 2 p of the de- 
composition order n max , the scale size /?, the flux F, the compo- 
nents of the ellipticity e, and the axis ratio r of the true shapelet 
model (denoted by 0) and the reconstruction using the two im- 
plementations. (R s ) is the mean distance in shapelet space be- 
tween the true model and the reconstructions (defined in Eq. 
([Toll). In the last line the mean centroid shift between the model 
and the reconstruction is given; consider that the IDL imple- 
mentation iteratively corrects the centroid, whereas the C++ 
implementation keeps it fixed. 



Model C++ 



IDL 



(F) 

^(e°,6 2 ) 
rl(r°,r) 

(R s ) 



1 

8.478 
1 
1 
1 
1 





9.791 
0.984 
8.478 
0.999 
0.999 
0.999 
0.989 
0.025 



11.969 

0.949 

8.406 

0.811 

0.894 

0.976 

0.771 

0.040 



00 



5 40 

s-r 20 



% 40 
20 



C + + 



IDL 



0.0 



0.5 



1.0 

AR r = \Ax r 



1.5 



2.0 



<|Ax c |> 







0.580 0.549 



Figure 14. Distribution of centroid errors for simulated images 
with underlying shapelet models. The centroid is obtained from 
the optimized shapelet coefficients. Only decompositions with 

x 2 



+ < 1 wAx 2 ' ,dl < 1 are considered. 



the expected impact of pixelization. It shows that the proce- 
dure of fixing the centroid from the pixel data (as suggested 
in sect. 14.4b produces errors which are not significantly larger 
than when the centroid is an optimized parameter, even when 
the galaxy model is a shapelet model. 

Thus, the underestimation of F in the IDL results must 
have different reasons. The distribution of the relative error 
6/3 = A/3/(3 {) = (J3 -jS )//? in Fig. [TJ shows again the larger 
scatter for /?"" but also reveals its slight underestimation, e.g. 
the maximum bin is at 6f3' dl = -0.01, which could give rise to 
an underestimated F. 

The distribution of 6f3 in the C++ results has roughly a 
Gaussian shape in the center but broad wings at both sides. 
When we set the <rpj intervals such that they correspond to 
68%, 95% and 99% confidence limits, we obtain cr^ ++ — 
(0.02,0.03,0.10) and af = (0.05,0.29,0.38). cr^ = 0.10 is 
the value we employed for the estimation of A/3 in sect. 14.61 

Since we know the coefficients of the underlying shapelet 
model 7°, we can quantify the error of the obtained shapelet 
coefficients in shapelet space by the simple Euclidean metric 



^ = £(7«-/ n ) 2 
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Figure 15. Distribution of relative ft errors for simulated im- 
ages with underlying shapelet models. Only decompositions 
with x 2 ' c++ ^ 1 an d X 2 '' dl ^ 1 are considered. The vertical 
lines indicate the 68% (solid, red), 95% (dotted, blue) and 99% 
(dashed, green) confidence limits. For the IDL results, the first 
and the last bin contains undershoot and overshoot, respec- 
tively. 



where we regard any uncertainty in the external parameters 
as error of the coefficients. Two important conclusions can 
be drawn from the corresponding numbers listed in Table Q] 
The goodness-of-fit in real space does not tell much about the 
goodness-of-fit in shapelet space (as all galaxies considered 
have^- 2 < 1); and the dominant source of coefficient error is 
not the pixel noise, otherwise the mean distance (R s ) should 
be close to the noise rms times the number of coefficients, 
(R s ) ~ 0.009. The realistic error estimate (which takes the addi- 
tional uncertainty in /3 into account, cf . sect. 14. 6t can reproduce 
such values of R*. 



5.4. Performance 

We ran the code on different machines. Our C++ implementa- 
tion ran on an ordinary desktop machine (Intel® Pentium® 4, 3 
GHz), whereas the IDL implementation ran on a more powerful 
machine (Intel® Xeon™, 2.8 GHz). To see how the platform 
impacts the runtime, we performed the easy-to-use BYTEmark 
benchmarlQ- a collection of realistic integer, floating point and 
memory problems - which indicates that the Xeon™ machine 
has an performance gain of factor 1.15 for floating point opera- 



7 http : // www . byte . com/bmark/bmark . htm 
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Table 2. Runtime comparison of the implementations: C++ on 
Pentium® 4 (3 GHz), IDL on Xeon™(2.8 GHz). The last col- 
umn shows a rough estimate of the runtime of the IDL imple- 
mentation on the slower Pentium® machine using the conver- 
sion factor 1.15. 
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« 607 min 
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ss 764 min 



tions; memory access and integer operations are equivalent on 
both machines. 

After calibrating the strength of the different machines, 
it was interesting to see what performance benefit could be 
gained by the use of C++ (see Tabled. 

Although our code tests more strictly whether the set of 
shapelet parameters used for the decomposition actually obeys 
Eq. ( fTBT l and thus needs more iterations, it outperforms the IDL 
code significantly, showing speedup factors between 5 and 14. 
Most of the performance gain is probably due to the fact that, in 
terms of computation time, IDL is not well suited for moderate 
to large numerical problems. 

6. Conclusions 

The shapelet technique provides a very powerful tool to de- 
scribe intrinsically smooth, compact objects. It incorporates the 
effects of noise such that it is able to extract all significant phys- 
ical information from the object. Quantities like the flux or the 
ellipticity can then be computed efficiently in the much reduced 
shapelet space. 

However, the trust in the estimators of the physical proper- 
ties relies on the assumption that there is a single decomposi- 
tion result. Since the shapelet decomposition depends on four 
external parameters (the scale size /3, the maximum decompo- 
sition order n max and the components of the centroid x c ), it is 
inevitable to choose appropriate parameters. Unfortunately, it 
is not sufficient for finding a reliable parameter set that the re- 
construction has residuals at noise level, i.e. x 2 ~ L In fact, 
increasing n max makes it more and more likely to find a con- 
tinuous and growing range of parameters obeying^- 2 < 1. This 
degeneracy has to be broken by ensuring that both conditions 
(fTTT i and ( fT3T l are fulfilled. If this is not the case, deviations of 
the shapelet parameters from the optimal values will introduce 
arbitrary errors on the shapelet coefficients and therefore on all 
quantities derived from them, even though^- 2 < 1. 

In addition, the sampled shapelet basis vectors have to re- 
main orthonormal to guarantee the correctness of the decompo- 
sition results. The orthonormality can be violated by undersam- 



pling, boundary truncation and introduction of non-orthogonal 
elements. There is no way to avoid undersampling, but the 
use of suitably large frames and the subtraction of noise be- 
fore the decomposition can remedy the other reasons for non- 
orthonormality. 

We showed that the proposed solutions do indeed result in 
a faithful representation of the decomposed objects, not only 
at the visual level but also concerning their physical properties. 
By comparing the results of our C++ implementation with the 
publicly available IDL code, we found that our code is both 
more reliable and has considerable performance benefits. 
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